import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from meteostat import Point, Daily
from datetime import datetime, timedelta
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller, kpss
import warnings
import seaborn as sns
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split, RandomizedSearchCV
from sklearn.metrics import mean_absolute_error, mean_squared_error
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.seasonal import STL
from sklearn.model_selection import GridSearchCV
from prophet import Prophet
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from scipy.stats import genextreme as gev, norm
from scipy.stats import genpareto

TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
Importing plotly failed. Interactive plots will not work.

# Zeitraum für Vergangenheitsdaten (letzte 50 Jahre, falls verfügbar)
start_date = datetime(1980, 1, 1)  # Statt 1940, um Lücken zu vermeiden
end_date = datetime(2024, 3, 1)

# Punkt für Bakersfield, USA
location = Point(35.3733, -119.0187)

# Daten von Meteostat API abrufen
data = Daily(location, start_date, end_date).fetch()

# Koordinaten von Bakersfield
lat, lon = 35.3733, -119.0187

# Karte von ganz Nordamerika
fig = plt.figure(figsize=(12, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([-170, -50, 5, 75], crs=ccrs.PlateCarree())  # Nordamerika-Ausschnitt

# Basis-Features
ax.add_feature(cfeature.LAND, facecolor='lightgray')
ax.add_feature(cfeature.OCEAN, facecolor='lightblue')
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle='-', linewidth=0.5)
ax.add_feature(cfeature.LAKES, alpha=0.5)
ax.add_feature(cfeature.RIVERS)

# Bakersfield markieren
ax.plot(lon, lat, marker='o', color='red', markersize=7, transform=ccrs.PlateCarree())
ax.text(lon + 2, lat, 'Bakersfield', transform=ccrs.PlateCarree(), fontsize=12)

# Titel
plt.title('Lage von Bakersfield in Nordamerika', fontsize=14)
plt.tight_layout()
plt.show()

print(data.describe())

print('-------------------------------------------------------------')

print(data.isnull().sum())

print('-------------------------------------------------------------')

data.hist(figsize=(15, 10))
plt.show()

print('-------------------------------------------------------------')

plt.figure(figsize=(10, 8))
correlation_matrix = data.corr()
sns.heatmap(correlation_matrix, annot=True)
plt.show()

print('-------------------------------------------------------------')

plt.figure(figsize=(10, 4))
sns.boxplot(x=data['tavg'])
plt.title('Boxplot Temperatur (tavg)')
plt.xlabel('Temperatur (°C)')
plt.show()

Q1 = data['tavg'].quantile(0.25)
Q3 = data['tavg'].quantile(0.75)
IQR = Q3 - Q1
lower_bound = Q1 - 1.5 * IQR
upper_bound = Q3 + 1.5 * IQR
data['tavg'] = data['tavg'].clip(lower=lower_bound, upper=upper_bound)

               tavg          tmin          tmax          prcp          snow  \
count  15808.000000  16132.000000  16132.000000  16132.000000  16130.000000   
mean      18.692371     11.926141     25.665485      0.435148      0.004712   
std        7.924751      6.987263      8.914693      2.099032      0.598407   
min      -17.800000     -7.200000      1.100000      0.000000      0.000000   
25%       12.200000      6.100000     18.300000      0.000000      0.000000   
50%       18.000000     11.100000     25.600000      0.000000      0.000000   
75%       25.600000     17.800000     33.300000      0.000000      0.000000   
max       37.200000     30.000000     46.100000     38.900000     76.000000   

              wdir          wspd         wpgt          pres    tsun  
count  6303.000000  16114.000000  5765.000000  15691.000000  1137.0  
mean    236.082818      9.953736    27.918508   1015.376573     0.0  
std     130.245126      3.673261     9.454404      5.222059     0.0  
min       0.000000      0.000000    11.200000    986.500000     0.0  
25%      94.000000      7.600000    22.300000   1011.500000     0.0  
50%     318.000000      9.700000    25.900000   1014.600000     0.0  
75%     333.000000     11.900000    31.300000   1019.000000     0.0  
max     360.000000     46.100000    92.500000   1034.000000     0.0  
-------------------------------------------------------------
tavg      324
tmin        0
tmax        0
prcp        0
snow        2
wdir     9829
wspd       18
wpgt    10367
pres      441
tsun    14995
dtype: int64
-------------------------------------------------------------

-------------------------------------------------------------

-------------------------------------------------------------

plt.figure(figsize=(15, 5))
plt.plot(data.index, data['tavg'], label='Average Temperature')
plt.xlabel('Date')
plt.ylabel('Temperature (°C)')
plt.title('Average Temperature over time')
plt.legend()
plt.show()

data_clean['tavg'].interpolate(method='time')

# 1. compute a calendar‐month z-score
monthly_mean = data['tavg'].groupby(data.index.month).transform('mean')
monthly_std  = data['tavg'].groupby(data.index.month).transform('std')
z_scores     = (data['tavg'] - monthly_mean) / monthly_std

# 2. flag >3σ outliers
threshold   = 3.0
is_outlier  = z_scores.abs() > threshold
outlier_dates = data.index[is_outlier]

# 3. group consecutive outlier‐dates into runs
bad_periods = []
if len(outlier_dates) > 0:
    run_start = outlier_dates[0]
    prev_date = outlier_dates[0]
    for curr_date in outlier_dates[1:]:
        if curr_date - prev_date > timedelta(days=1):
            # gap >1 day → end of run
            bad_periods.append((run_start, prev_date))
            run_start = curr_date
        prev_date = curr_date
    # final run
    bad_periods.append((run_start, prev_date))

print("Detected bad_periods:", bad_periods)

# 4. build a mask: outlier AND within one of those runs
mask = pd.Series(False, index=data.index)
idx  = data.index.to_series()
for start, end in bad_periods:
    mask |= is_outlier & idx.between(start, end)

# 5. mask as NaN
data_clean = data.copy()
data_clean.loc[mask, 'tavg'] = np.nan

# 6a. impute by time‐based linear interpolation
data_clean['tavg_interp'] = data_clean['tavg'].interpolate(method='time')

# 6b. alternative: rolling 7‐day median
data_clean['tavg_rollmed'] = (
    data_clean['tavg']
    .fillna(data_clean['tavg']
    .rolling('7D', center=True)
    .median()
    )
)

# 7. compare via plot
plt.figure(figsize=(15,5))
plt.plot(data['tavg'],         color='lightgray', label='raw')
plt.plot(data_clean['tavg_interp'],   label='interp (time)')
plt.plot(data_clean['tavg_rollmed'],  label='rolling median (7D)')
plt.legend()
plt.title("Raw vs. Imputed Average Temperature")
plt.xlabel("Date")
plt.ylabel("Temperature (°C)")
plt.show()

Detected bad_periods: [(Timestamp('1981-01-22 00:00:00'), Timestamp('1981-01-22 00:00:00')), (Timestamp('1989-02-05 00:00:00'), Timestamp('1989-02-06 00:00:00')), (Timestamp('1990-12-22 00:00:00'), Timestamp('1990-12-22 00:00:00')), (Timestamp('1990-12-24 00:00:00'), Timestamp('1990-12-24 00:00:00')), (Timestamp('1998-12-22 00:00:00'), Timestamp('1998-12-22 00:00:00')), (Timestamp('2001-01-23 00:00:00'), Timestamp('2001-01-23 00:00:00')), (Timestamp('2001-09-24 00:00:00'), Timestamp('2001-09-24 00:00:00')), (Timestamp('2002-10-31 00:00:00'), Timestamp('2002-11-01 00:00:00')), (Timestamp('2003-07-20 00:00:00'), Timestamp('2003-07-20 00:00:00')), (Timestamp('2003-08-31 00:00:00'), Timestamp('2003-08-31 00:00:00')), (Timestamp('2003-09-14 00:00:00'), Timestamp('2003-09-14 00:00:00')), (Timestamp('2003-12-02 00:00:00'), Timestamp('2003-12-02 00:00:00')), (Timestamp('2004-03-20 00:00:00'), Timestamp('2004-03-20 00:00:00')), (Timestamp('2010-11-05 00:00:00'), Timestamp('2010-11-05 00:00:00')), (Timestamp('2014-01-24 00:00:00'), Timestamp('2014-01-24 00:00:00')), (Timestamp('2017-01-28 00:00:00'), Timestamp('2017-01-28 00:00:00')), (Timestamp('2020-10-26 00:00:00'), Timestamp('2020-10-26 00:00:00'))]

data['tavg'] = data_clean['tavg_rollmed']

warnings.filterwarnings("ignore")
data.index = pd.to_datetime(data.index)
data = data.sort_index()
data = data.dropna(subset=['tavg'])  # STL benötigt vollständige Daten

# -----------------------------------------------------------------------------
# 2. STL mit optimierten LOESS-Fenstern
# -----------------------------------------------------------------------------
# period=365    → Jahreszyklus
# seasonal=365  → ganzjährige Saisonalität glätten (odd integer ≥ period)
# trend=1095    → 3-Jahres-Fenster für langfristigen Klimatrend
# low_pass=365  → Filter zur Entkopplung von Trend und Saisonalität
# robust=True   → resistent gegen Ausreisser im Fit
stl = STL(
    data['tavg'],
    period=365,
    seasonal=365,
    trend=1095,
    low_pass=367,
    robust=True
)
decomp = stl.fit()

# -----------------------------------------------------------------------------
# 3. Visualisierung
# -----------------------------------------------------------------------------
fig, axes = plt.subplots(4, 1, figsize=(12, 10), sharex=True)

# Original
axes[0].plot(data.index, data['tavg'], color='tab:blue')
axes[0].set_title('Originale Zeitreihe')
axes[0].set_ylabel('°C')

# Trend
axes[1].plot(data.index, decomp.trend, color='tab:red')
axes[1].set_title('Trendkomponente (LOESS-Fenster = 3 Jahre)')
axes[1].set_ylabel('°C')

# Saisonal
axes[2].plot(data.index, decomp.seasonal, color='tab:green')
axes[2].set_title('Saisonale Komponente (LOESS-Fenster = 1 Jahr)')
axes[2].set_ylabel('°C')

# Residuen
axes[3].scatter(data.index, decomp.resid, color='gray', alpha=0.4, s=5)
axes[3].set_title('Residuen (robust gefittet)')
axes[3].set_ylabel('°C')
axes[3].set_xlabel('Datum')

plt.tight_layout()
plt.show()

# 1) STL-Zerlegung mit optimierten LOESS-Fenstern
# -----------------------------------------------------------------------------
#   period=365    → Jahreszyklus
#   seasonal=365  → 1-Jahres-Fenster (ungerade ≥ period)
#   trend=1095    → 3-Jahres-Fenster für einen sehr glatten Klimatrend
#   low_pass=367  → ungerade & > period, trennt Trend↔Saison sauber
#   robust=True   → resistent gegen Ausreisser
stl = STL(
    data['tavg'],
    period=365,
    seasonal=365,
    trend=1095,
    low_pass=367,
    robust=True
)
decomp = stl.fit()

# -----------------------------------------------------------------------------
# 2) Nur die Trend-Komponente plotten
# -----------------------------------------------------------------------------
plt.figure(figsize=(12, 5))
plt.plot(
    decomp.trend.index,
    decomp.trend,
    color='tab:red',
    linewidth=2
)
plt.title('Optimierter Klimatrend der Temperatur (STL-Trend)')
plt.xlabel('Datum')
plt.ylabel('Temperatur (°C)')
plt.grid(True, linestyle='--', alpha=0.6)
plt.tight_layout()
plt.show()

# Fehlende Werte behandeln
data['tavg'] = data['tavg'].interpolate(method='linear')

window = 30  # Fensterbreite in Tagen
rolling_mean = data['tavg'].rolling(window).mean()
rolling_var  = data['tavg'].rolling(window).var()

plt.figure(figsize=(10,4))
plt.plot(data['tavg'], label='Original')
plt.plot(rolling_mean, label=f'{window}-Tage rollendes Mittel')
plt.fill_between(
    data.index,
    rolling_mean - np.sqrt(rolling_var),
    rolling_mean + np.sqrt(rolling_var),
    color='lightgray',
    alpha=0.5,
    label='±1 StdDev rollierend'
)
plt.title('Rolling Mean & Variance')
plt.legend()
plt.show()

# Statistische Tests auf Stationarität
# -----------------------------------------------------------------------------
# Augmented Dickey-Fuller (ADF)
adf_result = adfuller(data['tavg'], autolag='AIC')
print(f"ADF-Test:    p-Wert = {adf_result[1]:.4f}")

# KPSS
kpss_result = kpss(data['tavg'], nlags="auto", regression='c')
print(f"KPSS-Test:   p-Wert = {kpss_result[1]:.4f}")

# Entscheidung & ggf. Transformation
if (adf_result[1] < 0.05) and (kpss_result[1] >= 0.05):
    print("Zeitreihe gilt als stationär (ADF + KPSS erfüllt).")
    stationary = data['tavg']
else:
    print("Nicht stationär → Differenzierung der Zeitreihe.")
    stationary = data['tavg'].diff().dropna()

# Autokorrelationsanalyse (ACF & PACF)
fig, axes = plt.subplots(1, 2, figsize=(14,4))

plot_acf(stationary, lags=50, ax=axes[0])
axes[0].set_title("Autokorrelationsfunktion (ACF)")

plot_pacf(stationary, lags=50, ax=axes[1], method='ywm')
axes[1].set_title("Partielle Autokorrelationsfunktion (PACF)")

plt.tight_layout()
plt.show()

ADF-Test:    p-Wert = 0.0000
KPSS-Test:   p-Wert = 0.1000
Zeitreihe gilt als stationär (ADF + KPSS erfüllt).

# 2) Augmented Dickey-Fuller Test
adf_res = adfuller(data['tavg'].dropna(), autolag='AIC')
print("ADF-Statistik:", adf_res[0])
print("p-Wert:", adf_res[1])
for key, val in adf_res[4].items():
    print(f"Critial Value ({key}): {val}")

# 3) KPSS-Test (Nullhypothese: Stationär)
kpss_res = kpss(data['tavg'].dropna(), regression='c', nlags='auto')
print("\nKPSS-Statistik:", kpss_res[0])
print("p-Wert:", kpss_res[1])
for key, val in kpss_res[3].items():
    print(f"Critial Value ({key}): {val}")

ADF-Statistik: -10.051987486903496
p-Wert: 1.4046662561202967e-17
Critial Value (1%): -3.430761007049236
Critial Value (5%): -2.8617216481009167
Critial Value (10%): -2.5668666867669416

KPSS-Statistik: 0.13871588035233817
p-Wert: 0.1
Critial Value (10%): 0.347
Critial Value (5%): 0.463
Critial Value (2.5%): 0.574
Critial Value (1%): 0.739

# STL-Zerlegung zum Entfernen von Trend & Saison
stl = STL(
    data['tavg'],
    period=365,
    seasonal=365,
    trend=1095,
    low_pass=367,
    robust=True
)
res = stl.fit()
stationary_series = res.resid.dropna()

# Stationaritätstest der Residuen
adf_p  = adfuller(stationary_series)[1]
kpss_p = kpss(stationary_series, regression='c', nlags='auto')[1]
print(f"Residuen: ADF-p = {adf_p:.4f}, KPSS-p = {kpss_p:.4f}")

# Jahresmaxima der residualisierten Reihe
ann_max = stationary_series.resample('A').max().dropna()

# GEV-Anpassung
#    Achtung: scipy.stats.genextreme nutzt c = –ξ
c, loc, scale = gev.fit(ann_max)
print("Geschätzte GEV-Parameter (auf Residuen):")
print(f"  Shape (c = –ξ): {c:.4f}")
print(f"  Location (μ):   {loc:.2f}")
print(f"  Scale (σ):      {scale:.2f}")

# QQ-Plot
emp_q = np.sort(ann_max.values)
prob = (np.arange(1, len(emp_q)+1) - 0.5) / len(emp_q)
theo_q = gev.ppf(prob, c, loc=loc, scale=scale)

plt.figure(figsize=(6,4))
plt.scatter(theo_q, emp_q, alpha=0.7)
plt.plot([emp_q.min(), emp_q.max()],
         [emp_q.min(), emp_q.max()], 'k--')
plt.xlabel('Theoretische Quantile (GEV)')
plt.ylabel('Empirische Jahresmaxima')
plt.title('QQ-Plot GEV vs. Residual-Maxima')
plt.grid(True)
plt.show()

# Return-Level-Plot
return_periods = np.array([2,5,10,20,50,100])
prob_exc       = 1 - 1/return_periods
rl = gev.ppf(prob_exc, c, loc=loc, scale=scale)

# Rückübersetzung: Residual → Originalwert
# Da keine Zeitstempel für die Return Levels existieren, nehmen wir als typischen Wert den Mittelwert von Trend und Saisonalität der STL-Zerlegung.
# So kann eine typische Rückübersetzung auf das ursprüngliche Level erfolgen.

# Typischer Trend- und Saisonalitätswert (über alle Zeitpunkte gemittelt)
mean_trend = res.trend.mean()
mean_seasonal = res.seasonal.mean()
typical_level = mean_trend + mean_seasonal

# Rückübersetzte Return Levels
rl_orig = rl + typical_level

plt.figure(figsize=(6,4))
plt.plot(return_periods, rl, marker='o', label='Residuum')
plt.plot(return_periods, rl_orig, marker='s', label='Rückübersetzt (typisch)')
plt.xscale('log')
plt.xlabel('Rückkehrperiode (Jahre)')
plt.ylabel('Temperatur (°C)')
plt.title('Return-Level Plot: Residuen-GEV & Rückübersetzung')
plt.grid(which='both', ls='--')
plt.legend()
plt.show()

Residuen: ADF-p = 0.0000, KPSS-p = 0.1000
Geschätzte GEV-Parameter (auf Residuen):
  Shape (c = –ξ): 0.0676
  Location (μ):   10.32
  Scale (σ):      2.34

# Annual-Maxima aus den STL-Residuen (ann_max aus vorherigem Block)
# 1) Parametrischer Bootstrap
np.random.seed(42)
n_boot = 500
return_periods = np.array([2,5,10,20,50,100])
probs = 1 - 1/return_periods

params_bs = np.empty((n_boot, 3))
rl_bs     = np.empty((n_boot, len(return_periods)))

for i in range(n_boot):
    # simulierte Residual-Maxima
    sim = gev.rvs(c, loc=loc, scale=scale, size=len(ann_max))
    # Fit
    c_bs, loc_bs, scale_bs = gev.fit(sim)
    params_bs[i] = [c_bs, loc_bs, scale_bs]
    # Return Levels simuliert
    rl_bs[i] = gev.ppf(probs, c_bs, loc=loc_bs, scale=scale_bs)

# 95%-Konfidenzintervalle
ci_params = np.percentile(params_bs, [2.5,97.5], axis=0)
ci_rl     = np.percentile(rl_bs, [2.5,97.5], axis=0)

print("95%-CI für GEV-Parameter [shape, loc, scale]:")
print(ci_params)
print("\n95%-CI für Return Levels:")
for T, (lo, hi) in zip(return_periods, ci_rl.T):
    print(f" {T}-Jahre: {lo:.2f} – {hi:.2f}")

# 2) Mann–Kendall-Test auf ann_max
def mann_kendall(x):
    n = len(x)
    s = sum(np.sign(x[j] - x[i]) for i in range(n) for j in range(i+1, n))
    # Varianz mit Bindungen
    unique, counts = np.unique(x, return_counts=True)
    tie = counts[counts>1]
    var_s = (n*(n-1)*(2*n+5) - sum(t*(t-1)*(2*t+5) for t in tie)) / 18
    z = (s - np.sign(s)) / np.sqrt(var_s) if s != 0 else 0
    p = 2 * (1 - norm.cdf(abs(z)))
    return s, z, p

S, Z, p_val = mann_kendall(ann_max.values)
print(f"\nMann–Kendall auf Jahresmaxima: S={S}, Z={Z:.2f}, p={p_val:.3f}")

# 3) Declustering auf residualisierten tmax über 95%-Schwelle
threshold = stationary_series.quantile(0.95)
exceed = stationary_series > threshold
# Cluster-IDs generieren
cluster_id = (exceed & ~exceed.shift(1).fillna(False)).cumsum()
cluster_max = stationary_series[exceed].groupby(cluster_id).max()

print(f"\nAnzahl Cluster-Maxima: {len(cluster_max)}")
c_d, loc_d, scale_d = gev.fit(cluster_max)
print("GEV nach Declustering:")
print(f" c={c_d:.4f}, loc={loc_d:.2f}, scale={scale_d:.2f}")

# Histogramm-Vergleich
plt.figure(figsize=(6,4))
plt.hist(ann_max,        bins=20, alpha=0.6, label='Ann. Residual-Maxima')
plt.hist(cluster_max,    bins=20, alpha=0.6, label='Declustered-Maxima')
plt.legend()
plt.title('Jahres- vs. Cluster-Maxima')
plt.xlabel('Residual-Wert (°C)')
plt.ylabel('Anzahl')
plt.show()

95%-CI für GEV-Parameter [shape, loc, scale]:
[[-0.11918667  9.61385344  1.74449256]
 [ 0.33812847 11.20802015  2.88443463]]

95%-CI für Return Levels:
 2-Jahre: 10.41 – 12.09
 5-Jahre: 12.53 – 14.73
 10-Jahre: 13.76 – 16.43
 20-Jahre: 14.60 – 18.49
 50-Jahre: 15.58 – 21.38
 100-Jahre: 16.05 – 23.79

Mann–Kendall auf Jahresmaxima: S=-166.0, Z=-1.61, p=0.107

Anzahl Cluster-Maxima: 236
GEV nach Declustering:
 c=-0.6371, loc=9.30, scale=0.93

# für den Beispiel-Run hier:
np.random.seed(0)
dates = pd.date_range("1980-01-01", "2022-12-31", freq="D")
t = pd.Series(15 + 10*np.sin(2*np.pi*dates.dayofyear/365) 
              + np.random.randn(len(dates)), index=dates, name="tavg")


# 1) GEV-Fit auf jährliche Maxima & Risiko-Kennzahlen
# Jahresmaxima extrahieren
ann_max = t.resample("A").max()

# GEV anpassen (c = –ξ in scipy)
c_gev, loc_gev, scale_gev = gev.fit(ann_max)
print(f"GEV:    ξ = {-c_gev:.4f}, μ = {loc_gev:.2f}°C, σ = {scale_gev:.2f}°C")

# ── a) Wahrscheinlichkeit, dass JAHRESMAX > 40°C
thr = 40
p_exceed = 1 - gev.cdf(thr, c_gev, loc=loc_gev, scale=scale_gev)
print(f"P(Jahresmax > {thr}°C) = {p_exceed:.3%}")

# ── b) Expected Shortfall E[X | X > thr]  (für ξ < 1)
if c_gev < 1:
    # E[X|X>u] = u + (σ + ξ (u−μ)) / (1−ξ)
    es_gev = thr + (scale_gev + c_gev*(thr - loc_gev)) / (1 - c_gev)
    print(f"ES[X|X>{thr}°C] = {es_gev:.2f}°C")

# ── c) Simulation künftiger Jahresmaxima
sim_years = gev.rvs(c_gev, loc=loc_gev, scale=scale_gev, size=1000)
plt.figure(figsize=(6,4))
plt.hist(sim_years, bins=30, density=True, alpha=0.6, label="sim. Jahresmaxima")
x = np.linspace(sim_years.min(), sim_years.max(), 200)
plt.plot(x, gev.pdf(x, c_gev, loc=loc_gev, scale=scale_gev), 'k-', label="GEV-PDF")
plt.axvline(thr, color='r', ls='--', label=f"Schwelle {thr}°C")
plt.legend()
plt.title("Simulation von Jahresmaxima & GEV-Dichte")
plt.show()


# 2) POT/GPD auf tägliche Exzesse & Finanz-Risk-Kennzahlen
# Schwellenwahl: 95%-Quantil
u = t.quantile(0.95)
exc = t[t > u] - u

# GPD anpassen (loc=0 fixiert)
c_gpd, loc_gpd, scale_gpd = genpareto.fit(exc, floc=0)
print(f"GPD:    ξ = {c_gpd:.4f}, σ = {scale_gpd:.2f}°C")

# ── a) VaR und Expected Shortfall auf 99%-Niveau
q = 0.99
var_pot = u + genpareto.ppf(q, c_gpd, loc=0, scale=scale_gpd)
# ES potenziell: ES = u + σ/ (1−ξ) * (1−q)^(−ξ) − 1)/(1−ξ)
if c_gpd < 1:
    es_pot = u + (scale_gpd/(1 - c_gpd)) * ((1 - q)**(-c_gpd) - 1)
else:
    es_pot = np.nan
print(f"POT-VaR 99% = {var_pot:.2f}°C, POT-ES 99% = {es_pot:.2f}°C")

# ── b) QQ-Plot der Exzesse gegen GPD
emp = np.sort(exc.values)
prob = (np.arange(1,len(emp)+1) - 0.5) / len(emp)
theo = genpareto.ppf(prob, c_gpd, loc=0, scale=scale_gpd)

plt.figure(figsize=(6,4))
plt.scatter(theo, emp, alpha=0.7)
m, M = emp.min(), emp.max()
plt.plot([m,M],[m,M],'k--')
plt.xlabel("Theoretische Quantile (GPD)")
plt.ylabel("Empirische Exzesse")
plt.title("QQ-Plot GPD")
plt.grid(True)
plt.show()

# ── c) POT Return-Level-Plot (Jahre)
years = (t.index.max()-t.index.min()).days/365.25
λ = len(exc)/years
rp = np.array([2,5,10,20,50,100])
rl_pot = u + (scale_gpd/c_gpd)*((λ*rp)**c_gpd - 1)

plt.figure(figsize=(6,4))
plt.plot(rp, rl_pot, 's-')
plt.xscale('log')
plt.xlabel("Rückkehrperiode (Jahre)")
plt.ylabel("Rückkehrniveau (°C)")
plt.title("POT/GPD Return-Level-Plot")
plt.grid(which='both', ls='--')
plt.show()

GEV:    ξ = 0.0561, μ = 26.79°C, σ = 0.36°C
P(Jahresmax > 40°C) = 0.000%
ES[X|X>40°C] = 39.63°C

GPD:    ξ = -0.2114, σ = 0.81°C
POT-VaR 99% = 27.52°C, POT-ES 99% = 24.71°C

t = data['tavg']

# 1) Mean Residual Life Plot (Schwellenwahl) 
qs = np.linspace(0.85, 0.99, 15)
u_vals = t.quantile(qs)
mean_exc = [(t[t>u] - u).mean() for u in u_vals]

plt.figure(figsize=(6,4))
plt.plot(u_vals, mean_exc, 'o-')
plt.xlabel('Schwellenwert $u$ (°C)')
plt.ylabel('Mean Excess $E[X-u|X>u]$')
plt.title('Mean Residual Life Plot')
plt.grid(True)
plt.show()

# 2) Threshold‐Stability Plot 
xis, sigmas = [], []
for u in u_vals:
    exc = t[t>u] - u
    c, loc, scale = genpareto.fit(exc, floc=0)
    xis.append(c)
    sigmas.append(scale)

plt.figure(figsize=(6,4))
plt.plot(u_vals, xis, 'o-')
plt.xlabel('Schwellenwert $u$ (°C)')
plt.ylabel('geschätztes Shape ξ')
plt.title('Threshold Stability Plot')
plt.grid(True)
plt.show()

# 3) Exzesse extrahieren & Declustering 
u0 = t.quantile(0.95)     # wähle hier deinen Schwellenwert
exceed = t > u0
# Jeder neue Cluster beginnt, wenn exceed von False auf True wechselt
cluster_id = (exceed & ~exceed.shift(1).fillna(False)).cumsum()
# Cluster‐Maxima über dem Schwellenwert minus u0 = die Exzesse
exceedances = t[exceed].groupby(cluster_id).max() - u0

# 4) GPD-Fit auf die (declusterten) Exzesse 
c, loc, scale = genpareto.fit(exceedances, floc=0)
print(f"GPD-Parameter:  ξ={c:.4f},  scale={scale:.2f}  (loc={loc})")

# 5) QQ-Plot der Exzesse gegen GPD 
emp = np.sort(exceedances.values)
prob = (np.arange(1,len(emp)+1) - 0.5) / len(emp)
theo = genpareto.ppf(prob, c, loc=0, scale=scale)

plt.figure(figsize=(6,4))
plt.scatter(theo, emp)
m = min(emp.min(), theo.min()); M = max(emp.max(), theo.max())
plt.plot([m,M],[m,M],'k--')
plt.xlabel('Theoretische Quantile (GPD)')
plt.ylabel('Empirische Exzesse')
plt.title('QQ-Plot GPD')
plt.grid(True)
plt.show()

# 6) Return-Level-Plot (Jahres-Rückkehrperioden) 
# Jahresrate an Exzessen:
years = (data.index.max() - data.index.min()).days / 365.25
lambda_year = len(exceedances) / years

return_periods = np.array([2,5,10,20,50,100])  # in Jahren
# Formel: x_T = u0 + (σ/ξ)[ (λ_year·T)^ξ − 1 ]
rl = u0 + (scale / c) * ((lambda_year * return_periods)**c - 1)

plt.figure(figsize=(6,4))
plt.plot(return_periods, rl, 'o-')
plt.xscale('log')
plt.xlabel('Rückkehrperiode (Jahre)')
plt.ylabel('Rückkehrniveau (°C)')
plt.title('Return-Level-Plot (POT/GPD)')
plt.grid(True, which='both', ls='--')
plt.show()

GPD-Parameter:  ξ=-0.3500,  scale=2.34  (loc=0)

# Nur Residuen ohne NaNs (wegen Trendberechnung bei STL)
resid = decomp.resid.dropna()

# Trainings-/Testsplit
train_size = int(len(resid) * 0.8)
resid_train = resid.iloc[:train_size]
resid_test = resid.iloc[train_size:]

# SARIMA-Modell auf die Residuen trainieren (monatlich saisonal, Periode 12)
sarima_resid_model = SARIMAX(resid_train, 
                             order=(1, 0, 1), 
                             seasonal_order=(1, 0, 1, 12),
                             enforce_stationarity=False, 
                             enforce_invertibility=False)
sarima_resid_result = sarima_resid_model.fit(disp=False)

# Residuen-Vorhersage für Testzeitraum
resid_forecast = sarima_resid_result.forecast(steps=len(resid_test))
resid_forecast.index = resid_test.index  # Index angleichen

# Trend und Saisonalität passend zum Testzeitraum (nur wo resid_test gültig ist)
trend = decomp.trend.loc[resid_test.index]
seasonal = decomp.seasonal.loc[resid_test.index]

# Rekonstruktion = Trend + Saisonalität + prognostizierte Residuen
reconstructed_forecast = trend + seasonal + resid_forecast

# Tatsächliche Werte zum Vergleich
actual = data['tavg'].loc[reconstructed_forecast.index]

# Plot
plt.figure(figsize=(14, 5))
plt.plot(actual.index, actual, label="Original", color='black', alpha=0.6)
plt.plot(reconstructed_forecast.index, reconstructed_forecast, label="Rekonstruktion (SARIMA auf Residuen)", color='darkred')
plt.xlabel("Datum")
plt.ylabel("Temperatur (°C)")
plt.title("SARIMA-Vorhersage (monatlich saisonal) auf STL-Residuen")
plt.legend()
plt.grid(True, linestyle='--', alpha=0.6)
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()

# MAE ausgeben
mae_resid = mean_absolute_error(actual, reconstructed_forecast)
print(f"MAE der rekonstruierten Vorhersage: {mae_resid:.2f} °C")

MAE der rekonstruierten Vorhersage: 3.34 °C

data['tavg'] = data_clean['tavg_rollmed']

# Datum als Index setzen
data.index = pd.to_datetime(data.index)

# Fehlende Werte interpolieren
# data['tavg'] = data['tavg'].interpolate(method='time')

# Prophet erwartet zwei Spalten: ds (Datum), y (Wert)
df_prophet = data[['tavg']].reset_index()
df_prophet.columns = ['ds', 'y']

# Prophet-Modell aufsetzen
model = Prophet(
    yearly_seasonality=True,
    daily_seasonality=False,
    weekly_seasonality=False,
    seasonality_mode='additive'
)
model.fit(df_prophet)

# Zukunftsdaten erzeugen (10 Jahre)
future = model.make_future_dataframe(periods=10*365)
forecast = model.predict(future)

# Farbenblind-freundliche Farben (Okabe & Ito Palette)
colors = {
    'historical': '#0072B2',  # Blau
    'forecast': '#E69F00',    # Orange
    'uncertainty': '#999999'  # Grau
}

# Plotten
fig, ax = plt.subplots(figsize=(14, 6))

# Historische Temperatur
ax.plot(df_prophet['ds'], df_prophet['y'], label="Vergangene Temperatur",
        color=colors['historical'], linestyle='-', linewidth=2)

# Prognose
ax.plot(forecast['ds'], forecast['yhat'], label="Prophet Vorhersage",
        color=colors['forecast'], linestyle='--', linewidth=2)

# Unsicherheitsbereich
ax.fill_between(forecast['ds'], forecast['yhat_lower'], forecast['yhat_upper'],
                color=colors['uncertainty'], alpha=0.3, label="Unsicherheitsbereich")

# Achsenbeschriftungen und Titel
ax.set_xlabel("Datum")
ax.set_ylabel("Temperatur (°C)")
ax.set_title("Langfristige Temperaturvorhersage für Bakersfield (Prophet-Modell)")

# Legende, Gitter und Layout
ax.legend()
ax.grid(True, linestyle='--', alpha=0.6)
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()

12:34:07 - cmdstanpy - INFO - Chain [1] start processing
12:34:08 - cmdstanpy - INFO - Chain [1] done processing

# Datum als Index setzen
data.index = pd.to_datetime(data.index)

# Fehlende Werte interpolieren
# data['tavg'] = data['tavg'].interpolate(method='linear')

# Prophet erwartet zwei Spalten: ds (Datum), y (Wert)
df_prophet = data[['tavg']].reset_index()
df_prophet.columns = ['ds', 'y']

# Prophet-Modell aufsetzen
model = Prophet(
    yearly_seasonality=True,
    daily_seasonality=False,
    weekly_seasonality=False,
    seasonality_mode='additive'
)
model.fit(df_prophet)

# Zukunftsdaten erzeugen (10 Jahre)
future = model.make_future_dataframe(periods=50*365)
forecast = model.predict(future)

# Farbenblind-freundliche Farben (Okabe & Ito Palette)
colors = {
    'historical': '#0072B2',  # Blau
    'forecast': '#E69F00',    # Orange
    'uncertainty': '#999999'  # Grau
}

# Plotten
fig, ax = plt.subplots(figsize=(14, 6))

# Historische Temperatur
ax.plot(df_prophet['ds'], df_prophet['y'], label="Vergangene Temperatur",
        color=colors['historical'], linestyle='-', linewidth=2)

# Prognose
ax.plot(forecast['ds'], forecast['yhat'], label="Prophet Vorhersage",
        color=colors['forecast'], linestyle='--', linewidth=2)

# Unsicherheitsbereich
ax.fill_between(forecast['ds'], forecast['yhat_lower'], forecast['yhat_upper'],
                color=colors['uncertainty'], alpha=0.3, label="Unsicherheitsbereich")

# Achsenbeschriftungen und Titel
ax.set_xlabel("Datum")
ax.set_ylabel("Temperatur (°C)")
ax.set_title("Langfristige Temperaturvorhersage für Bakersfield (Prophet-Modell)")

# Legende, Gitter und Layout
ax.legend()
ax.grid(True, linestyle='--', alpha=0.6)
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()

12:34:11 - cmdstanpy - INFO - Chain [1] start processing
12:34:12 - cmdstanpy - INFO - Chain [1] done processing

warnings.filterwarnings("ignore")

# Sicherstellen, dass das Datum der Index ist
data.index = pd.to_datetime(data.index)

# Überprüfen, ob Meteostat-Daten verfügbar sind
print("Verfügbarer Zeitraum:", data.index.min(), "bis", data.index.max())

# Fehlende Werte in der Temperaturspalte behandeln
# data['tavg'] = data['tavg'].interpolate(method='linear')  # Interpolieren
# data['tavg'] = data['tavg'].fillna(method='bfill').fillna(method='ffill')  # Backfill + Forwardfill als zusätzliche Sicherung

# Prüfen, ob noch NaN-Werte vorhanden sind
print("Fehlende Werte nach Verarbeitung:", data['tavg'].isna().sum())

# Fehlende Werte in weiteren Wettervariablen auffüllen
data[['tmin', 'tmax', 'prcp', 'wspd', 'pres']] = data[['tmin', 'tmax', 'prcp', 'wspd', 'pres']].fillna(method='ffill').fillna(method='bfill')

# Temperatur-Zeitreihe extrahieren
temperature_series = data['tavg'].dropna()

# Trainings- und Testsplit
train_size = int(len(temperature_series) * 0.8)
train, test = temperature_series.iloc[:train_size], temperature_series.iloc[train_size:]

# SARIMA-Modell trainieren
sarima_model = SARIMAX(train, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12), enforce_stationarity=False, enforce_invertibility=False)
sarima_result = sarima_model.fit()

# Langfristige Vorhersage für 5 Jahre
forecast_steps = 5 * 365  # 5 Jahre
future_dates = pd.date_range(start=end_date, periods=forecast_steps, freq='D')
sarima_forecast = sarima_result.forecast(steps=forecast_steps)

# Random Forest Modell mit zusätzlichen Wetterfaktoren
features = data[['tmin', 'tmax', 'prcp', 'wspd', 'pres']].fillna(method='ffill').fillna(method='bfill')
target = data['tavg'].fillna(method='ffill').fillna(method='bfill')

# Überprüfung der Datenverteilung im Trainingsdatensatz
train_stats = target.describe()

# Trainings- und Testdaten für Random Forest
X_train, X_test, y_train, y_test = train_test_split(features, target, test_size=0.2, shuffle=False)

# Random Forest Modell optimieren
rf_model = RandomForestRegressor(n_estimators=500, max_depth=20, random_state=42)
rf_model.fit(X_train, y_train)

# Überprüfung: Sind future_features konstant?
future_features = features.rolling(window=30, min_periods=1).mean().iloc[-forecast_steps:].reset_index(drop=True)

# Falls zu wenige Werte vorhanden sind, mit Durchschnittswerten auffüllen
if future_features.shape[0] < forecast_steps:
    additional_rows = pd.DataFrame(np.tile(future_features.mean().values, (forecast_steps - future_features.shape[0], 1)), columns=future_features.columns)
    future_features = pd.concat([future_features, additional_rows], ignore_index=True)

# Random Forest Vorhersage
rf_forecast = rf_model.predict(future_features)

# Berechnung von Fehlerwerten
sarima_mae = mean_absolute_error(test, sarima_result.forecast(steps=len(test)))
rf_mae = mean_absolute_error(y_test, rf_model.predict(X_test))

# Ergebnisse visualisieren
plt.figure(figsize=(14, 6))
plt.plot(temperature_series.index, temperature_series, label="Vergangene Temperatur", color="black", alpha=0.6)
plt.plot(future_dates, sarima_forecast, label=f"SARIMA Vorhersage (MAE={sarima_mae:.2f})", color="blue", linestyle="--", linewidth=2)
plt.plot(future_dates, rf_forecast, label=f"Random Forest Vorhersage (MAE={rf_mae:.2f})", color="orange", linestyle="-.", linewidth=2)
plt.xlabel("Datum")
plt.ylabel("Temperatur (°C)")
plt.title("Langfristige Temperaturvorhersage für Bakersfield (SARIMA & Random Forest)")
plt.legend()
plt.grid(True, linestyle='--', alpha=0.6)
plt.xticks(rotation=45)
plt.show()

# Rückgabe der wichtigsten Werte
train_stats, sarima_mae, rf_mae

Verfügbarer Zeitraum: 1980-01-01 00:00:00 bis 2024-03-01 00:00:00
Fehlende Werte nach Verarbeitung: 0

(count    15957.000000
 mean        18.738247
 std          7.879740
 min         -0.600000
 25%         12.200000
 50%         18.000000
 75%         25.700000
 max         37.200000
 Name: tavg, dtype: float64,
 7.7163296875724425,
 0.6983166408206647)

warnings.filterwarnings("ignore")

# Sicherstellen, dass das Datum der Index ist
data.index = pd.to_datetime(data.index)

# Überprüfen, ob Meteostat-Daten verfügbar sind
print("Verfügbarer Zeitraum:", data.index.min(), "bis", data.index.max())

# Fehlende Werte in der Temperaturspalte behandeln
# data['tavg'] = data['tavg'].interpolate(method='linear')  # Interpolieren
# data['tavg'] = data['tavg'].fillna(method='bfill').fillna(method='ffill')  # Backfill + Forwardfill als zusätzliche Sicherung

# Prüfen, ob noch NaN-Werte vorhanden sind
print("Fehlende Werte nach Verarbeitung:", data['tavg'].isna().sum())

# Fehlende Werte in weiteren Wettervariablen auffüllen
data[['tmin', 'tmax', 'prcp', 'wspd', 'pres']] = data[['tmin', 'tmax', 'prcp', 'wspd', 'pres']].fillna(method='ffill').fillna(method='bfill')

# Temperatur-Zeitreihe extrahieren
temperature_series = data['tavg'].dropna()

# Trainings- und Testsplit
train_size = int(len(temperature_series) * 0.8)
train, test = temperature_series.iloc[:train_size], temperature_series.iloc[train_size:]

# SARIMA-Modell trainieren
sarima_model = SARIMAX(train, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12), enforce_stationarity=False, enforce_invertibility=False)
sarima_result = sarima_model.fit()

# Langfristige Vorhersage für 15 Jahre
forecast_steps = 15 * 365
future_dates = pd.date_range(start=end_date, periods=forecast_steps, freq='D')
sarima_forecast = sarima_result.forecast(steps=forecast_steps)

# Random Forest Modell mit zusätzlichen Wetterfaktoren
features = data[['tmin', 'tmax', 'prcp', 'wspd', 'pres']].fillna(method='ffill').fillna(method='bfill')
target = data['tavg'].fillna(method='ffill').fillna(method='bfill')

# Überprüfung der Datenverteilung im Trainingsdatensatz
train_stats = target.describe()

# Trainings- und Testdaten für Random Forest
X_train, X_test, y_train, y_test = train_test_split(features, target, test_size=0.2, shuffle=False)

# Random Forest Modell optimieren
rf_model = RandomForestRegressor(n_estimators=500, max_depth=20, random_state=42)
rf_model.fit(X_train, y_train)

# Überprüfung: Sind future_features konstant?
future_features = features.rolling(window=30, min_periods=1).mean().iloc[-forecast_steps:].reset_index(drop=True)

# Falls zu wenige Werte vorhanden sind, mit Durchschnittswerten auffüllen
if future_features.shape[0] < forecast_steps:
    additional_rows = pd.DataFrame(np.tile(future_features.mean().values, (forecast_steps - future_features.shape[0], 1)), columns=future_features.columns)
    future_features = pd.concat([future_features, additional_rows], ignore_index=True)

# Random Forest Vorhersage
rf_forecast = rf_model.predict(future_features)

# Berechnung von Fehlerwerten
sarima_mae = mean_absolute_error(test, sarima_result.forecast(steps=len(test)))
rf_mae = mean_absolute_error(y_test, rf_model.predict(X_test))

# Ergebnisse visualisieren
plt.figure(figsize=(14, 6))
plt.plot(temperature_series.index, temperature_series, label="Vergangene Temperatur", color="black", alpha=0.6)
plt.plot(future_dates, sarima_forecast, label=f"SARIMA Vorhersage (MAE={sarima_mae:.2f})", color="blue", linestyle="--", linewidth=2)
plt.plot(future_dates, rf_forecast, label=f"Random Forest Vorhersage (MAE={rf_mae:.2f})", color="orange", linestyle="-.", linewidth=2)
plt.xlabel("Datum")
plt.ylabel("Temperatur (°C)")
plt.title("Langfristige Temperaturvorhersage für Bakersfield (SARIMA & Random Forest)")
plt.legend()
plt.grid(True, linestyle='--', alpha=0.6)
plt.xticks(rotation=45)
plt.show()

# Rückgabe der wichtigsten Werte
train_stats, sarima_mae, rf_mae

Verfügbarer Zeitraum: 1980-01-01 00:00:00 bis 2024-03-01 00:00:00
Fehlende Werte nach Verarbeitung: 0

(count    15957.000000
 mean        18.738247
 std          7.879740
 min         -0.600000
 25%         12.200000
 50%         18.000000
 75%         25.700000
 max         37.200000
 Name: tavg, dtype: float64,
 7.7163296875724425,
 0.6983166408206647)

hybrid_forecast = sarima_forecast + rf_model.predict(future_features)

warnings.filterwarnings("ignore")

# Datum als Index sicherstellen
data.index = pd.to_datetime(data.index)

# Fehlende Werte behandeln
data['tavg'] = data['tavg'].fillna(method='bfill').fillna(method='ffill')
data[['tmin', 'tmax', 'prcp', 'wspd', 'pres']] = data[['tmin', 'tmax', 'prcp', 'wspd', 'pres']].fillna(method='ffill').fillna(method='bfill')

# Temperatur-Zeitreihe extrahieren
temperature_series = data['tavg'].dropna()

# Trainings- und Testsplit (80/20)
train_size = int(len(temperature_series) * 0.8)
train = temperature_series.iloc[:train_size]
test = temperature_series.iloc[train_size:]

# SARIMA trainieren
sarima_model = SARIMAX(train, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12),
                       enforce_stationarity=False, enforce_invertibility=False)
sarima_result = sarima_model.fit()

# Residuen berechnen
sarima_fitted = sarima_result.fittedvalues
residuals = train - sarima_fitted

# Input-Features für RF vorbereiten
features = data[['tmin', 'tmax', 'prcp', 'wspd', 'pres']]
features = features.fillna(method='ffill').fillna(method='bfill')

# RF auf Residuen trainieren (nur für SARIMA-Zeitraum)
X_rf = features.iloc[:len(residuals)]
y_rf = residuals

rf_model = RandomForestRegressor(n_estimators=500, max_depth=20, random_state=42)
rf_model.fit(X_rf, y_rf)

# Zukunft vorbereiten (5 Jahre)
forecast_steps = 5 * 365
end_date = data.index.max()
future_dates = pd.date_range(start=end_date + pd.Timedelta(days=1), periods=forecast_steps, freq='D')
sarima_forecast = sarima_result.forecast(steps=forecast_steps)

# Zukünftige Input-Features für RF (rolling average zur Glättung)
future_features = features.rolling(window=30, min_periods=1).mean().iloc[-forecast_steps:].reset_index(drop=True)

# Falls zu kurz: mit Mittelwerten auffüllen
if future_features.shape[0] < forecast_steps:
    filler = pd.DataFrame(np.tile(future_features.mean().values,
                                  (forecast_steps - future_features.shape[0], 1)),
                          columns=future_features.columns)
    future_features = pd.concat([future_features, filler], ignore_index=True)

# RF-Residuen-Vorhersage + Hybrid-Prognose
rf_residual_forecast = rf_model.predict(future_features)
hybrid_forecast = sarima_forecast + rf_residual_forecast

# MAE für SARIMA auf Testdaten
from sklearn.metrics import mean_absolute_error
sarima_test_forecast = sarima_result.forecast(steps=len(test))
sarima_mae = mean_absolute_error(test, sarima_test_forecast)

# Visualisierung
plt.figure(figsize=(14, 6))
plt.plot(temperature_series.index, temperature_series, label="Vergangene Temperatur", color="black", alpha=0.6)
plt.plot(future_dates, sarima_forecast, label=f"SARIMA Vorhersage (MAE={sarima_mae:.2f})", color="blue", linestyle="--", linewidth=2)
plt.plot(future_dates, hybrid_forecast, label="Hybride Vorhersage (SARIMA + RF)", color="green", linestyle="-.", linewidth=2)
plt.xlabel("Datum")
plt.ylabel("Temperatur (°C)")
plt.title("Langfristige Temperaturvorhersage für Bakersfield (Hybrides Modell)")
plt.legend()
plt.grid(True, linestyle='--', alpha=0.6)
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()

warnings.filterwarnings("ignore")

# Datum als Index sicherstellen
data.index = pd.to_datetime(data.index)

# Fehlende Werte behandeln
data['tavg'] = data['tavg'].fillna(method='bfill').fillna(method='ffill')
data[['tmin', 'tmax', 'prcp', 'wspd', 'pres']] = data[['tmin', 'tmax', 'prcp', 'wspd', 'pres']].fillna(method='ffill').fillna(method='bfill')

# Temperatur-Zeitreihe extrahieren
temperature_series = data['tavg'].dropna()

# Trainings- und Testsplit (80/20)
train_size = int(len(temperature_series) * 0.8)
train = temperature_series.iloc[:train_size]
test = temperature_series.iloc[train_size:]

# SARIMA trainieren
sarima_model = SARIMAX(train, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12),
                       enforce_stationarity=False, enforce_invertibility=False)
sarima_result = sarima_model.fit()

# Residuen berechnen
sarima_fitted = sarima_result.fittedvalues
residuals = train - sarima_fitted

# Input-Features für RF vorbereiten
features = data[['tmin', 'tmax', 'prcp', 'wspd', 'pres']]
features = features.fillna(method='ffill').fillna(method='bfill')

# RF auf Residuen trainieren (nur für SARIMA-Zeitraum)
X_rf = features.iloc[:len(residuals)]
y_rf = residuals

rf_model = RandomForestRegressor(n_estimators=500, max_depth=20, random_state=42)
rf_model.fit(X_rf, y_rf)

# Zukunft vorbereiten (5 Jahre)
forecast_steps = 15 * 365
end_date = data.index.max()
future_dates = pd.date_range(start=end_date + pd.Timedelta(days=1), periods=forecast_steps, freq='D')
sarima_forecast = sarima_result.forecast(steps=forecast_steps)

# Zukünftige Input-Features für RF (rolling average zur Glättung)
future_features = features.rolling(window=30, min_periods=1).mean().iloc[-forecast_steps:].reset_index(drop=True)

# Falls zu kurz: mit Mittelwerten auffüllen
if future_features.shape[0] < forecast_steps:
    filler = pd.DataFrame(np.tile(future_features.mean().values,
                                  (forecast_steps - future_features.shape[0], 1)),
                          columns=future_features.columns)
    future_features = pd.concat([future_features, filler], ignore_index=True)

# RF-Residuen-Vorhersage + Hybrid-Prognose
rf_residual_forecast = rf_model.predict(future_features)
hybrid_forecast = sarima_forecast + rf_residual_forecast

# MAE für SARIMA auf Testdaten
sarima_test_forecast = sarima_result.forecast(steps=len(test))
sarima_mae = mean_absolute_error(test, sarima_test_forecast)

# Visualisierung
plt.figure(figsize=(14, 6))
plt.plot(temperature_series.index, temperature_series, label="Vergangene Temperatur", color="black", alpha=0.6)
plt.plot(future_dates, sarima_forecast, label=f"SARIMA Vorhersage (MAE={sarima_mae:.2f})", color="blue", linestyle="--", linewidth=2)
plt.plot(future_dates, hybrid_forecast, label="Hybride Vorhersage (SARIMA + RF)", color="green", linestyle="-.", linewidth=2)
plt.xlabel("Datum")
plt.ylabel("Temperatur (°C)")
plt.title("Langfristige Temperaturvorhersage für Bakersfield (Hybrides Modell)")
plt.legend()
plt.grid(True, linestyle='--', alpha=0.6)
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()

# 1) SARIMA-Vorhersage auf Testdaten
sarima_test_forecast = sarima_result.forecast(steps=len(test))

# 2) Hybrid-Vorhersage auf Testdaten
#    – SARIMA-Forecast plus RF-Residuen-Forecast derselben Länge
sarima_part = sarima_result.forecast(steps=len(test))
#    Für die RF-Residuen braucht es Input-Features für denselben Zeitraum:
features_test = features.iloc[train_size:train_size+len(test)]
rf_resid_test = rf_model.predict(features_test)
hybrid_test_forecast = sarima_part + rf_resid_test

# 3) Kennzahlen berechnen
def mape(y_true, y_pred):
    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100

# SARIMA
sarima_mae  = mean_absolute_error(test, sarima_test_forecast)
sarima_rmse = np.sqrt(mean_squared_error(test, sarima_test_forecast))
sarima_mape = mape(test, sarima_test_forecast)

# Hybrid
hybrid_mae  = mean_absolute_error(test, hybrid_test_forecast)
hybrid_rmse = np.sqrt(mean_squared_error(test, hybrid_test_forecast))
hybrid_mape = mape(test, hybrid_test_forecast)

print("SARIMA       MAE: {:.3f}, RMSE: {:.3f}, MAPE: {:.2f}%".format(
    sarima_mae, sarima_rmse, sarima_mape))
print("Hybrid SAR+RF MAE: {:.3f}, RMSE: {:.3f}, MAPE: {:.2f}%".format(
    hybrid_mae, hybrid_rmse, hybrid_mape))

SARIMA       MAE: 7.716, RMSE: 8.966, MAPE: nan%
Hybrid SAR+RF MAE: 7.604, RMSE: 8.837, MAPE: nan%

plt.figure(figsize=(12,5))
plt.plot(test.index, test, label="Echte Temperatur", color="black")
plt.plot(test.index, sarima_test_forecast, label="SARIMA", linestyle="--")
plt.plot(test.index, hybrid_test_forecast, label="Hybrid SAR+RF", linestyle="-.")
plt.legend()
plt.xlabel("Datum")
plt.ylabel("Temperatur (°C)")
plt.title("Modellvergleich auf Testdaten")
plt.show()

Komponente	Aufgabe im hybriden Modell
SARIMA	Modelliert Trend + Grundsaisonalität auf Basis der Temperatur-Zeitreihe
Random Forest	Lernt aus Residuen des SARIMA-Modells + nutzt externe Einflussgrössen
Hybridmodell	Addiert die durch RF vorhergesagten Residuen zur SARIMA-Prognose

Kennzahl	SARIMA	Hybrid (SARIMA + RF)	Differenz (Hybrid–SARIMA)
MAE	7.716 °C	7.604 °C	–0.112 °C (besser)
RMSE	8.966 °C	8.837 °C	–0.129 °C (besser)
MAPE	nan %	nan %	–

SPZ - OMC¶

Abruf historischer Wetterdaten für Bakersfield (1980–2024)¶

Explorative Datenanalyse (EDA) der Wetterzeitreihe¶

Statistische Kennzahlen¶

Prüfung auf fehlende Werte¶

Histogramme¶

Korrelationsmatrix¶

Boxplot der Durchschnittstemperatur (tavg)¶

Bereinigung von Ausreissern mittels IQR-Methode¶

Zeitliche Entwicklung der durchschnittlichen Temperatur¶

Interpretation:¶

Umgang mit Ausreissern und fehlenden Werten: Zeitbasierte Imputation¶

1. Z-Score-basiertes Outlier-Detection-Verfahren (monatsspezifisch)¶

2. Identifikation zusammenhängender "Bad Periods"¶

3. Maskieren der Ausreisser in tavg¶

4. Imputation fehlender Werte¶

Saisonale Dekomposition der Temperaturzeitreihe mittels STL¶

1. Mathematische Grundlage der STL-Zerlegung¶

2. Implementierung mit STL (Seasonal-Trend Decomposition using Loess)¶

3. Grafische Darstellung der Komponenten¶

4. Interpretation & Anwendung¶

Extraktion und Visualisierung des Temperaturtrends mittels saisonaler Dekomposition¶

1. Mathematische Struktur der additiven Dekomposition¶

2. Implementierung mit seasonal_decompose¶

3. Visualisierung der Trendkomponente¶

4. Anwendung¶

Stationaritätsprüfung und Korrelationsanalyse zur Modellvorbereitung (ADF-Test, ACF, PACF)¶

1. Vorbereitung der Zeitreihe¶

2. Stationaritätsprüfung mit Augmented-Dickey-Fuller-Test (ADF)¶

3. Analyse der Autokorrelationen (ACF und PACF)¶

ACF – Autokorrelationsfunktion:¶

PACF – Partielle Autokorrelationsfunktion:¶

4. Fazit¶

1. Überblick der Stationaritätsprüfung¶

2. ADF-Test¶

3. KPSS-Test¶

4. Visuelle Diagnose (ACF/PACF)¶

5. Kritische Diskussion¶

Extremwert Analyse¶

Erster Ansatz¶

Block-Maxima-Methode¶

Auswertung der Extremwert-Analyse¶

1. Unsicherheitsabschätzung (Parametrischer Bootstrap)¶

2. Mann–Kendall-Trendtest auf Jahresmaxima¶

3. Declustering von Hitze-Ereignissen¶

4. Vergleich: Jahresmaxima vs. Declustered Maxima¶

5. Empfehlungen und nächste Schritte¶

Peaks-Over-Thershold (POT)-Methode¶

Wozu können Extremwertverteilungen genutzt werden?¶

1. Klimaforschung & Umweltanalyse¶

2. Technische Infrastruktur & Katastrophenschutz¶

3. Finanz- und Versicherungswesen¶

4. Qualitätskontrolle & technische Prozesse¶

5. Ökosysteme & Umwelttoxikologie¶

6. Simulation & Unsicherheitsanalyse¶

Fazit¶

Interpretation der POT-Analyse (Peaks-Over-Threshold)¶

1. Mean Residual Life Plot¶

2. Threshold Stability Plot¶

3. QQ-Plot der Exzesse gegen GPD¶

4. Return-Level-Plot (POT/GPD)¶

SARIMA-Modellierung der STL-Residuen (monatliche Saisonalität) mit Rekonstruktion der Temperatur¶

1. Modellidee: Hybrider Ansatz¶

2. Modellkonfiguration¶

3. Rekonstruktion der Zeitreihe¶

4. Visualisierung & Modellbewertung¶

5. Fazit¶

Temperaturprognose mit Prophet-Modell (10 Jahre) — Mathematische Modellierung¶

1. Datenformatierung¶

2. Mathematische Struktur des Prophet-Modells¶

3. Saisonalitätskomponente (Fourierreihe)¶

4. Prognose und Unsicherheit¶

5. Visualisierung¶

6. Vorteil des Prophet-Modells¶

Vergleich von SARIMA und Random Forest zur Temperaturprognose (mathematisches Modell + Analyse)¶

1. SARIMA-Modell: Mathematische Struktur¶

2. Random Forest Regressor: Mathematischer Hintergrund¶

3. Visualisierung und Interpretation¶

4. Fazit¶

Welchen Einfluss hätte ein hybrider Ansatz auf die enttäuschenden SARIMA-Vorhersagen?¶

Boxplot der Durchschnittstemperatur (`tavg`)¶

3. Maskieren der Ausreisser in `tavg`¶

2. Implementierung mit `STL` (Seasonal-Trend Decomposition using Loess)¶

2. Implementierung mit `seasonal_decompose`¶